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I.  Introduction 

In  the  consideration  of  propagation  of  electromagnetic  waves  in 
inhomogeneous  anisotropic  stratified  media, Maxwell's  equations  reduce  to  a 
system  of  foTir  first  order  linear  ordinary  differential  equations.   The 
treatment  of  these  equations  has  resulted  in  the  present  extension  of  the 
Peano-Baker  method.  However,  since  higher  order  systems  of  this  type  occur 
frequently  in  all  hranches  of  applied  mathematics  the  theory  is  developed 
for  a  system  of  •%!■  equations.  The  application  of  the  analysis  to  propaga- 
tion in  the  ionosphere  will  appear  in  a  subsequent  report. 

Consider  the  system  of  ordinary  different  led  eou-^tions. 

1)  ^  =  A(z)  U(z) 

in  which  U(z)  ia  an  n^rowed  column  vector  and  a(z)  is  a  square  matrix  function  of  z. 
When  the  matrix  A  is  constant,  the  solutions  of  (l)  are  linear  cochin 3t ions  of 
exponential  functions  with  constant  coefficients  (if  the  latent  roots  of  A  are 
distinct)  which  are  interpreted  as  damped,  growing  or  constant  oscillations  or 
wares,  or  else  represent  non-oscillatoiy  growth  or  decay.  These  solutions  pxe   so 
useful  in  applications  that  it  is  natural  to  seek  similar  solutions  when  A  is  not 
constant.  Many  attempts  to  ohtaln  such  solutions  have  heen  made  and  much  of  this 
work  has  often  been  xuinecessarily  repeated  partly  because  there  is  no  general 
treatment  of  the  system  (l)  from  this  point  of  view. 

The  nearest  approach  to  such  a  treatment  is  contained  in  the  WKB  method 
which  yields  an  approximate  solution  of  a  linear  second  order  equation  as  a  sum 
of  exponentials  with  variable  coefficients  and  exponents.  In  some  present 'it ions 
of  the  method  higher  order  approximations  are  also  obtained,  but  it  is  generally 
not  clear  whether  the  successive  approximations  converge  to  a  solution  of  the 
equation  or  not.   For  a  large  class  of  second  order  equations  involving  a  parameter 
in  a  particular  way,  Birkhoff  and  Langa:  have  settled  this  matter  by  showing  that 
the  WKB  approximation  is  the  leading  term  in  an  asyntptotic  expansion  of  the  solution 
with  respect  to  the  parameter. 

In  all  other  cases,  however,  when  the  equation  or  sya'^em   does  not  involve 
a  parameter  (or  does  but  not  in  the  appropriate  manner)  this  result  isn't  appli- 
cable and  the  significance  of  the  V/KB  a-;)proximation  remains  in  doubt.  But  recently 
several  workers  have  shoxm  that  the  WKB  approximation  for  the  reduced  one  dimen- 
sional wave  equation  can  be  obtained  as  the  first  term  in  the  application  of  the 
Picard  iteration  method  to  an  equivalent  system  of  integral  equations   .Thus  when 


the  Picard  method  yields  a  solution,  the  \/KB  approximation  is  the  first  term  In.  a 
convergent  series.  The  successive  terms  in  the  series  have  heen  interpreted  in 
terms  of  multiple  reflections. 

For  systems  or  eqTiations  of  order  greater  tlian  two,  there  is  no  general 
analogue  of  the  V/KB  method,  although  for  equations  containing  a  parameter  in  the 
appropriate  manner  Birkhoff b  theorem  assures  the  existence  of  an  asymptotic  ex- 
paxLsion  with  respect  to  the  parameter.  The  leading  terms  in  the  expansion  are 
exponentials  of  -che  desired  type,  and  thus  provide  a  generalization  of  the  WKB 
method  to  such  equations.  But  for  all  those  cases  in  which  the  system  does  not 
involve  a  parameter,  there  is  no  generalisation  available. 

In  this  article  we  present  i/hat  seems  to  "be  a  new  approach  to  the 
original  problem  of  obtaining  exponential-type  solutions  of  the  system  (l).  We 
begin  with  the  Peano-Baker  matrizant  representation  of  the  solution  of  (l), 
which  is  known  to  be  uniformly  convergent  in  a  specified  large  region,  althou^ 
it  converges  very  slowly  and  therefore  doesn't  exhibit  the  behavior  of  the 
solution.  By  applying  certain  transfonnntions  to  this  representation  we  obtain 
a  representation  of  the  desired  type — the  solution  is  a  sum  of  e:q)onential8 
with  varifble  coefficients.  This  representation,  which  is  of  course  still  uni- 
formly convergent  in  a  large  region,  ig  rapidly  convergent  in  many  casee  and  further- 
more exhibits  clearly  the  behavior  of  the  solution.  Thus  one  of  the  v/eaknesses 
of  the  matrizant  method,  its  slow  convergence,  is  removed. 

When  the  result  is  specialized  to  a  second  order  system  equivalent  to 
a  single  second  order  equation,  the  first  term  in  the  representation  is  identi- 
cal with  the  '.rt3  approximation.   Thus  this  representation  may  be  considered  to 
be  the  natural  generalization  of  the  WKB  method  to  higher  order  systems.   In 
addition  it  shows  that  the  'rfKB  approximation  is  the  first  term  of  a  uniformly 
convergent  representation  of  the  solution,  and  the  region  of  convergence  is 
known.  The  representation  fails,  and  must  be  modified,  at  points  where  the 
elementary  divisors  of  A(z)  are  not  simple  i.e.  at  some  points  where  two  or 
more  latent  roots  (eigenvalues)  of  a(z)  are  equal.   In  the  second  order  case 
these  are  the  ordinary  turninig-points,  and  thus  we  may  also  call  them  turning- 
points  in  the  higher  order  case.  The  original  matrizant  representation  is 
valid  in  a  neighborhood  of  each  such  point,  and  thus  provides  the  connection  be- 
tween the  other  Tepresentations  of  the  rigorouB  solution  in  its  vicinity. 


If  the  system  (l)  involves  a  parameter  in  an  appropriate  way,  then  we 
can  deduce  an  asymptotic  expansion  of  the  solution  with  respect  to  this  parameter, 
from  our  representation.   Thus  we  can  essentially  deduce  Birkhoff's  theorem 
concernins;  the  existence  of  an  asymptotic  expansion  and  further  ottain  explicit 
expressions  for  the  successive  terms.   In  particular  the  first  term  is  still 
the  (generalized)  WKB  approximation,  and  thus  we  prorve  that  this  approximation 
is  Eilso  the  leading  term  in  an  asymptotic  expansion  of  the  solution  provided 
that  the  system  contains  a  parameter  appropriately.  The  connection  formulae 
for  the  asymptotic  expansion  are  not  considered. 

The  exponential  terms  can  be  readily  interpreted  as  oscillations  or 
waves,  just  es  in  the  case  where  A  is  constant.  Each  exponent  is  the  integral 
of  one  latent  root  of  A  and  the  coefficient  of  each  exponential  is  variable,  in 
fact  it  is  a  series  of  ntultiple  integrals.  These  integrals  have  an  immediate 
interpretation  in  terms  of  continuous  multiple  reflection  caused  by  the  vari- 
ation of  A(z),  Thus  it  is  generally  not  possible  for  a  single  wave  to  occur 
alone,  i.e.  to  yield  a  solution  itself.   Instead  any  one  v/ave  gives  rise  to 
the  others  by  the  continuous  reflections  just  mentioned.   In  addition  "lumped" 
reflections  may  occur  from  regions  between  the  turning  points  where  the  separate 
waves  may  be  strongly  "coupled"  together.  The  overall  effect  of  these  reflections 
Is  described  by  the  scattering  matrix,  which  is  defined  and  expressed  in  terms 
of  our  representstlon  of  the  solution. 
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II  The  Feano-Baker  Method 

We  wish  to  examine  the  solution  of  the  system  1)  for  real  and  com- 
plex yalues  of  z. 


^^  ^  =  A(0  U(z)  ;  U(z  )  =  U 

dz  0     0 

In  the  complex  case  we  require  that: 

(a)  The  n  elements  a  . (z)  of  the  coefficient  matrix  A(z)  are 
single  valued  and  analytic  in  the  "star"  region  formed  hy 

non-intersecting  straight  lines  from  all  the  singularities 

and  "branch  points  of  each  of  the  a.  .(z)  to  Infinity.  The 


For  z  real  we  require  instead  that: 

(h)  The  a.  .(z)  are  single  valued,  hounded  and  integrable  in  an 
interval  containing  the  point  z  . 
Under  either  of  the  ahove  conditions  it  has  heen  proved  hy  various  methods 
that  the  solution  of  the  one  point  hoijindary  value  prohlem  stated  "by  equa- 
tion 1)  exists  and  is  unique.   In  the  present  paper  we  wish  to  exploit  the 
method  of  Peano[ll  and  Baker  [2,3i^]  . 

Following  Baker  [4]  we  define  the  square  matrix  of  n-rows,  which 

* 

he  has  called  the  matrizant  of  A(z),  hy 

z  z       x^ 

2)      X^UCx)"!  =  I  +  yA(x^)dx^  +  /a(x^)  y  ACx^)  dx^  dx^  +  

^0  2^  z„        z 

0  0  0 


*  One  is  led  to  the  same  series  hy  applying  the  Picard  iteration  scheme  to 
the  system  (1).   Or  considering  the  equivalent  system  of  integral  equa- 
tions 

U(z)  =  U^  +  /  A(x)  U(x)  dx  , 
z 

0 

the  Neumann  expansion   of   the  solution  yields    the  series    (2 ) .      Also  Tolterra's 
[13]    defi!'itlon   of   the    "left   integral"  of  a  matrix   is   precisely  its  matri^.ant, 
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Here  I  is  the  n  x  n  unit  matrix.      In  the  complex  case  all   integrals   in  (2) 
are  over  the  same  path  vdiich  does  not   intersect  amy  of  the  barriers  of  condition 
(a)  bat   is  otherwise  arbitrary.     The  series   for   the  elements  of  /i    \a(x)  ^con- 
verge uniformly  and  absolutely  [5],   end  thus   term  by  term  differentiation  is 
valid  and  yields 

z  s 

3)  fe     /1(a(x)]    =     A(z)    nlA(x))     . 

Hence  the  solution  of  the   systemd)  is 

z 
1^)  U(z)  =    n    ^  A(x)"j    U^ 

0 

This  series  converges  uniformly  and  absolutely  throughout  the  entire  star 
region  of  (a)  or  in  the  interval  of  (b),  which  regions  are  much  larger  than 
those  obtained  by  most  other  methods.  This  feature  is  compensated  for  by  the 
fact  that  the  convergence  is  usually  so  slow  that  the  solution  in  this  form  is 
not  practical  for  computation. 
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III.     Transformation  of  the  Solution 

In  order  to  transform  the  solution  (U)  to  a  more  rapidly  convergent 
form  we  will  employ  the  following  identities: 

15)  A     [a{x)]    =     P(z)    n    \   p-\x)  A(x)  P(x)  -  P"^x)  ^^Jp-\z^ 

0  0    *■ 

55         /  -J  Z 

-HL     \  A(x)}  =     exp     I        a(x)  dx 

0  o^'olo  0  -^ 

z 
4"      [X(x)]    .       A°    [xCx)] 


6) 


7) 


9) 


4^    \Mx)\      =       n\AU)]    n    ^A(x)] 


Identity  (5)  holds  for  any  non-sing ilar  matrix  P(z)  which  is  analytic,  in  the 
complex  case,  and  has  a  hounded  integrahle  derivative,  in  the  real  case.   Identity 
(6)  holds  for  any  matrix  such  that  A(z')  ACz")  =  A(z")  A(z«)  for  all  z'  and  z" 
on  the  path  of  integration.   In  particular  this  applies  to  all  diagonal  and 
constant  matrices.  Identities  (7),  (8)  and  (9)  hold  for  arhitraiy  matrices 
provided  the  natrizants  exist.   Identity  (6)  is  proved  in  appendix  I;  the  others 
appear  in  the  literature  [  1+, 5.6j  . 

Applying  identity  (5),  the  solution  given  by  (U-)  may  be  written 

10)  U(z)  =  P(z)  A    \  P-^X)  A(x)  P(x)  -  P-^X)  Mil]  P-1(,^)  U^  . 

0  ^  ■' 

If  P(z)  can  he  chosen  such  that  i  P'-^AP  -  P~^  ^  [  is  diagonal,  identity  (6) 

leads  to  a  solution  in  closed  form.  However,  there  is  no  general  procedure 

for  detemining  such  a  matrix  P  and  therefore  we  merely  attempt  to  diagonalize 
the  first  term,  P~  AP, 
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At  any  point  z   this  is  possible  if  and  only  if  the  elementary 
divisors  of  the  matrix  (A(z)  -  M)  are  simple,  where  X  is  a  constant  [7], 
Then  we  have 


11) 


P'-'Cz)  A(z)  Hz)   = 


where  the  X.(z)  are  the  roots  of  the  characteristic  equation  J^  73 


12) 


dat  (A(z)  -  M)  =  0 


We  assTine  that  there  are  regions  in  the  star  of  condition  (a)  or  the  interval 
of  condition  (h)  throughout  which  this  diagonal! zat ion  is  possible  and  we  re- 
strict our  attention  to  such  regions. 

We  now  introduce  the  abbreviation 


13) 


dz 


Then  using  equations  (ll)  and  (I3)  we  may  write 


Vii 


lu)   |p-W-^g] 


V-^as 


o 


o 


n  nn 


/?  y 


"21 


^nl 


n-1 ,  n 
,n-l-0 


Here  the  r  Xz)     are  the  elements  of  R(z).  We  use  the  equation  (lU)  in  the 

solution,  equation  (lO)jand  expand  the  resulting  matrlzant  with  the  aid  of 
identities  (6)  and  (7).  This  yields 

15)       U(z)  «=  P(z)  •  exp  J    L(x)dx  .   Q   fM(z^,Tj|  .  P'^(z^)  U^   . 

0 

In  this   expression  L(x)  is   the   diagonal  matrix  on  the  right   side  of   equation   (lU) 

and  has   elements  a.(i)  h  >v  +r    ..     M(z   ,z)  is   the  matrix  whose  elements  are 
J  J     Jj  o 


16) 


ij 


(z)  5  (1-  (f,  J     r,Az)  /° 


'ir     ^ij' 


The  representf'-tion  of  the  solution  given  by  equation  (I5)  is   the  deBired  rapidly 
convergent  form  which  exhibits  its  behavior. 

In  order  to  er.press   the  solution  more   clearly  we  introduce  the  two 
column  matrices  C(z   ,z)  and  B(z   )  by  the  definitions 


17)  OU.z)  5     A{m(z^,x)]  B(z   ) 

'  0'  zt       0        '         0 


Kz  )  =    P'-'Cz  )    u 
0  00 


Farther,  the  elements  of  C(z  , z)  may  be  written  as 


18 


)        V^o»*^=  J^^  As  ^^o'^^^^^o^   '  ^  =  1.2 ^   • 


Here  the  B  (z  )  are  the  elements  of  the  matrix  B(z^)  and  the  i  i-   (z  z) 

so  O  Ij   0 

are  the  elements  of  the  matrizant  in  equation  (I7).  These  are,  by  the  definitions 
(2)  and  (16) 

*  n    '        f^ 


TTsing  equations   (I7)  in  the  solution  (15)  yields 


20) 


U(z)  =  P(z)  •    exp    (         L(x)  dx  •   C(z^,z)       . 


Thus  a  particular  component  of  this  solution  may  be  written  ae 


21)      u.(z)  =  Yl 
^  J=l 


i    'JJ^^^^J   ^^^^^ 


P^j(z)  e 


Cj(z^,z)  :  1=1,2.  ....  n. 


Here  u.(z)  is  a  component  of  U(z).  the  ?.<(*)  are  the  elements  of  P(z),  and  the 
elements  of  L(x)  have  been  written  explicitly. 
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Finally,  by  taking  account  of  the  definition  of  ^^^C3c)  It  le 
T)OssnDle  to  write  the  solution  given  "by  equation  (21)  In  the  form 


22)      u.(z)  =  r 


/  qjj(x)dx 


Pi^^-^-  ° 


det  P(z) 


f   Xj(x) 


dx 


C.(z^,s)  •  det  P(z  )  , 
Jo'  0 


Here     q^^C^)  is   defined  "by 


23) 


<iQ,.(x) 


•I  n  do 

q,.(^)s     ^—      E      P..U)     -^ 


v:here  *jj_(x)  is  the  cofactor  of  P^jd)  in  det  P(x)  (see  appendix  Hi. 

Equation  22  yields  the  desired  representntion  of  the  solution.   Sa.ch 
component  of  the  solution  is  a  sum  of  exponential  functions  with  variable  ex- 
ponents and  variable  coefficients.  Each  exponent  is  the  integral  of  a  latent 
root  of  tlie  matrix  A.  Each  coefficient  is  a  product  of  a  function  in  closed 
fona  and  another,  C.,  which  is  an  infinite  series  of  maltiple  integrals.   If 
only  the  first  term,  a  constant,  in  each  G.  is  retained  the  solution  coincides 
with  the  VKB  approximation  in  the  second  order  case  for  vdiich  the  latter  is 
dsfined,  (see  section  7),  Thais  the  WKB  approximation  is  seen  to  be  the  first 
term  in  a  uniformly  and  absolutely  convergent  representation  of  the  solution, 
Jurthennore  the  leading  terms  in  equation  (22)  can  be  considered  a  natural 
generalization  of  the  WEB  approximation  to  higher  order  systems.  The  successive 
terras  in  C.  can  be  simply  interpreted  as  corresponding  to  multiple  continuous 
reflection  caused  by  the  variation  in  A  (see  section  V).   This  explains why  a 
single  exponential  term  (or  a  single  wave)  cannot  yield  a  solution  by  itself: 
each  wave  produces  other  waves  by  continuous  reflection,   '^nly  v/hen  the  deriva- 
tives of  A  are  zero  (i.e,  when  A  is  constant)  are  these  reflections  absent, 
and  thus  in  this  case  a  single  wave  can  constitute  a  solution. 
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IV.  Connection  Across  Turning  Points. 

The  representation  of  the  solution  given  Isy  equation  (22)  Is  not  valid 
if  there  are  points  on  the  path  of  integration  at  which  A  cannot  "be  dia^onalized, 
and  in  the  second  order  case  these  are  Just  the  turning  points  of  the  ordinary 
WKB  method.  Thus  we  will  also  call  them  turning  points  in  the  general  case,  and 
we  will  asBTiiue  that  they  are  isolated.   In  the  conplex  case  turning  points  offer 
no  difficulty,  because  it  is  only  necessarj'  to  restrict  the  path  of  integration 
to  avoid  them,  as  was  done  in  the  preceding  section.  Since  any  point  wloich  is 
not  Itself  a  turning  point  can  be  reached  from  the  initial  point  z     by  such  a 
path,  the  representation  given  by  equation  (22)  holds  at  all  points  except  turning 
points.  In  the  case  of  z  real,  however,  the  representation  In  equation  (22) 
holds  only  when  z  and  z  are  not  separated  by  a  turning  point.  To  extend  the 
solution  across  a  turning  point,  or  to  coxmect  the  representations  on  opposite 
sides  of  a  turning  point,  it  is  necessary  to  eimjloy  a  different  representation 
of  the  solution,   'Ve  will  use  the  original  matrlzant  representation  for  this 
purpose. 

Thus,  suppose  there  are  turning  points  only  in  the  interval  z^  <  z  <  Zj^, 
and  we  wish  to  obtain  the  solution  at  f  >  z^.  By  the  method  of  the  preceding 
section  this  may  be  v/rltten  as 

2h)  U(^)  =  P(<)  exp  j    L(y)dy  O  ^  M(z3^,y)]p-^(z^)  U(z^) 

^1  ^ 

where  all  matrices  are  defined  as  above.  However,  this  solution  involves  U(zj_) 
and  we  are  given  only  U( z  )  =  U^.  Applying  the  solution  in  the  form  given  by 
equation  (U)  we  obtain 

25)  U(zi)-  n     {A(y)|  U^   . 

Now  the  solution  for  ^  >z^   may  be  v;rltten  as^ 
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where  the  colTinm  matrix  D(z. ,  ^  )  is  defined  by 


27)  ^U,.^)S  n  |M(z^.y)]p-^(.^)  -Q^  |A(jr)j  U^ 

Lution  (2£ 
JAj(y)dy 


As  before  a  particular  component  of  the  solution  (26)  is 


2s)    u.(n  =  2I 


Pjj^^^-' 


det  P(  ^  ) 


Dj(zj^,^)«det  P(zj^)  ;  1=1,2 n 


Equations  (27)  and  (28 )  give  the  solution  for  f  >  2,  while  equations  (20)-(22) 
now  apply  only  for  z  <   z^.  The  connection  across  the  interval  containing  the 
turning  points  has  been  famished  by  the  matrizant  appearing  in  equation  (25). 
The  solution  for  ^  >  Zj_  may  be  obtained  from  that  for  z  <  2  by  determining  the 
square  matrix  T  (  ^  ,z, ;  z,z  )  such  that 

29)  DCzj^,  P  =  T  (^.Zj^;  z.z^)  C(z^,z)   . 
Prom  equations  (I7)  and  (27)  this  matrix  is  found  to  be 

30)  T(^.z^;z.z^)5  n  |M(z^.y)]  P-^(z^)  ij  |A(y)]p(z^)  -H^^MCz^.y)! 


and  this  furnishes  the  connection  between  any  points  z  and    ^   in  their  respective 
ranges. 
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V,   Scattering  Matrix 

The  connection  matrix  T  defined  by  equation  (30)  relates  the  amplitudes 
of  the  waves  at  any  two  points  ^  and  z,  when  the  phases  of  the  waves  are  re- 
ferred to  the  points  z,  and  z  respectively.   This  matrix  can,  of  course,  be  de- 
fined whether  or  not  turning  points  are  present;  if  there  are  none  it  becomes 
simpler  since  z  and  z,  may  then  be  identified.   It  is  frequently  useful  to 
define  another  matrix  related  to  T,  the  scattering  matrix  S.   In  order  to  de- 
fine S  it  is  first  necessary  to  designate  some  of  the  waves  at  f  and  some  at  z 
as  "incoming"  or  "incident"  waves;  the  remainder  are  called  "outgoing"  or 
"scattered"  waves.  The  scattering  matrix  yields  the  amplitudes  of  the  scattered 
waves  in  terms  of  those  of  the  incident  waves.  For  this  to  be  possible  it  is 
necessary  that  half  of  the  2n  waves  (n  each  at  |"  and  z)  be  Incoming  and  the 
other  half  outgoing. 

Let  us  suppose,  then,  that  the  roots  are  so  adered  that  the  first  k 


waves  at  z  and  the  last  n-k  waves  at  f  are  Incident  waves. 

C, 


we  may  write 


C(5 


)    = 


(U) 


.k+l 


31) 


D(z,.P 


D-^ 


Vi 


■■/ 


Here  C^  and  J)^^  are  the  amplitudes  of  the  incident  waves  jC"^  and  if  are  the 
aJBplitudes  of  the  scattered  waves.  Equation  (29)  may  now  be  written  in  block 
form  as 


32) 


-n-k 


^kk 


n-k.k 


\ 

^k.n-k        * 


n-k,  n-k 


.n-k 
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The  scattering  matrix  S  is  defined  "by 

33)  1  '   ' 

FroBi  the  last  two  equations  we  obtain 


3^^)  S(^,z^; 


n-k 


provided  that 

35)  -^-^Vk.n-k^^     • 

The  scattering  matrix  S  given  "by  eq,.  (3I+)  yields  the  amplitudes  of 
the  scattered  waves  at  ^  emd  z  in  terms  of  the  amplitudes  of  the  incident  waves 
at  the  same  points,  when  the  phases  of  the  wpves  at  f  and  z  are  referred  to 
z-^   and  Zq  respectively.   Its  definition  is  independent  of  the  presence  or  ahsence 
of  turning  points.   In  most  cases  of  interest  S  approaches  a  limit  as  ^  end  z 
approach  plus  and  minus  infinity  respectively,  and  this  limit  is  v/hat  is  usually 
called  the  scattering  matrix.   It  depends  only  upon  the  reference  points  z,  and 
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VI,  Asymptotic  Expansion  of  the  Solution. 

If  the  matrix  A(z)  depends  upon  a  parameter  k  in  an  appropriate  way, 
then  the  Bolution  of  equations  (l)  posseBses  an  asymptotic  expansion  in  reciprocal 
powers  of  k  in  regions  where  the  latent  roots  of  A  are  distinct.  This  result 
was  o'btained  by  Birklioff  [s  J  for  an  nth  order  equation  and  has  since  heen  ex- 
tended to  nth  order  systems.  The  connection  "between  the  asymptotic  expansions 
in  different  regions,  when  turning  points  are  present,  has  "been  treated  "by 
Langerfg"!  for  second  order  equations,  but  very  little  has  been  done  on  this 
problem  for  higher  order  systems.   From  the  viewpoint  of  the  present  paper,  the 
significance  of  Birkhoff 's  result  is  that  the  asymptotic  expansion  which  he  ob- 
tains is  just  a  sum  of  exponential  functions  with  variable  coefficients  and 
variable  exponents.  In  the  second  order  case  the  leading  terras  in  the  expansion 
coincide  with  the  WKB  approximation,  thus  proving  that  the  WKB  approximation 
yields  the  leading  term  in  an  asymptotic  expansion  of  the  solution  If  the 
equation  depends  upon  a  parameter  appropriately.  Similarly,  in  the  higher  order 
case  the  leading  terras  of  our  representation  of  the  solution  in  eq.  (?2)  coin- 
cide with  the  leading  terras  in  an  asymptotic  expansion  of  Birkhoff 's  type. 

Since  the  asymptotic  expansion  can  be  deduced  easily  from  our  form  of 
solution,  and  since  this  derivation  may  clarify  the  role  of  tlie  WKB  approxi- 
mation as  the  leading  term  in  both  convergent  and  asymptotic  expansions,  we  will 
derive  It.   Suppose,  then,  that  the  matrix  A  dapends  upon  a  parameter  k  such 
that  A(z,k)  =  k  A(z)  where  a(z)  is  Independent  of  k.  The  matrix  P(z)  is  then 
chosen  to  dlagonalize  A(z)  and  hence  it,  P~  (z)  and  R(z)  are  also  Independent 
of  k.  The  elements  of  the  matrix  M(z  ,z)  now  become,  from  equation  (l6), 

/  (Rjj(z)-a^^Cz)dzj  k|*  (X^(z)-X^(7j)dz 


36) 


.,j(k.  =  ).(l-<^,j) 


R,j(z) 


Here  the  X.Cz)  are  the  eigenvalues  of  the  matrix  A(z).  The  solution  of  the 


system 


dU(z) 
dz 


A(z,k)  U(z)  is  from  equation  (2l) 


37)    u,(z,k)  =  f2    P.  1^2^®  ° 
^       j=l  ^    '' 


f   H  (z)dz  1  kj  X.(z)dz 


Hji^^.^^^^^i*-" 


A.(^o. 


;k)B, 


i  =  1.2. 
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Tiie  nj_j(zQ,z;k)   are  o'btained  "by   replacing  the  '^^Az)  appearing  in  equation  (19) 
Ijy  the  mj^j(k,z)   of  equation   (36).      The  asymptotic   e.rpansion  of   this   solution 
is   determined  by  that  of  the  -Tl     .   (z   ,z;k)   . 
Let  U3   introduce   the   notation 


f,j(z)  .  (1-  5^j) 


/(Hj^j(z)-R^^(z))dz 


ij' 


38) 


z 
g     (z)  5/^     (Xj(z)-?v^(z))  dz 


with  §4  4^0^  ~  °  ^^"^  i  ?*  J«     Now  we  have 


39) 


n,j(z^.z;k)  =     ^^j  .f^   f^j(.^)   e'^^J^'^     dx, 
0 


r=l  ^  z 


We  wish  to  integrate  "by  parts  (m+l)  tiraee  in  eaich  of  the  above  terms; 
this  requires  that  the  f.  .(z)  and  g.  a(z)  "be  differ entiable  (■♦l)  and(m+2)  times 
respectively.   In  case  (a) this  condition  is  automatically  satisfied  but  case  (b) 

must  be  amended  to  include  differentiability  of  order  (m+2)  of  the  a.  .(z).  To 

ij 
prove  the  asymptotic  nature  of  the  resulting  expansion  we  must  restrict  our- 
selves to  that  portion  of  the  complex  plane  or  real  interval  where 

UO)  Ee  [kg^j(z)]  <  0 

for  all  g..(z).     But   from  the  definition  (3S)   it   is  seen  that  g^ 4( z)  =  -  gj^Cz) 

and  so  condition  (UO)  implies 


Kj(ol 


Under  the  above   conditions  we  may  now   integrate  by  parts   (m+l)   tines   in  the 
first  integral  of  equation  (39).     This  yields 
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0 

k   Jz       dy 


<^J\.,) 


where 


o<^^^u)=(^)^.(^)-^^^(^)=-(¥^) 


dg  (z) 
Here  we  have  used  the  relation  — ■^ =  X.(z)-X.  (z)  f   0,  which  is  otvious  from 

(38)  and  the  absence  of  turning  points.  By  condition  (UO)  it  is  now  apparent 

that  for  |k\  ->  00  the  first  integral  of  equation  (39)  has  the  asymptotic 

expansion  given  "by  the  first  ao  terms  on  the  right  of  equation  (^1),   Similarly 

it  Is  easily  demonstrated  that  the  p-th  iterated  integral  terms  in  (39)  have 

asymptotic  expansions  In  powers  of  7-  starting  v.dth  terms  of  —  .  By  adding 

all  of  these  expansions  we  finally  obtain  the  asymptotic  expansion  of 
_Q  .  .(z  ,2;k)  and  hence  of  the  solution  (37). 

The  terms  imiltiplying  the  second  square  bracket  of  equation  (37)  con- 
stitute what  we  have  called  in  the  general  case,  or  is  exactly  for  the  second 
order  eouation,  the  W.B.K,  solution.  Purther  these  square  brackets  have  asymptotic 
expansions  v;hich  consist  of  a  constant  plus  a  power  series  in  r-  .  Thus  it  1b 
shown  that  in  general  the  W.B.K.  solution  is  the  leading  term  in  the  asymptotic 
expansion  of  the  rigorous  solution. 


-17- 

VII,     Second  Order  Bcpiation. 

We   shall  now  apply  the  preceeding  analysis   to   the   second  order 
differential   equation 

dz 


with  the  one  point  houndary  conditions 


yl 


This  equation  with  the  houndary  conditions  may  "be  put  into  the  form  of  the 
system  (l)  hy  introducing  vi^iz)   «  '^'(z);  U2(z)  =  ^J^     •  and  the  matrices: 


m 


(0  1         \  /u,(z)\       /^ov 

A(z)-  I  I    U(z)  =         .  U  =/ 

V-^(0    2a(z)/  V^U)I         °     [^l/ 


To  solve  the  resulting  equation  by  the  method  of  sections  II   and  III  we  are 
led  to  consider  the  algehralc  equation 

XI  -  A(z)|  =  X^-2a(z)  >v+  tCz)  »=  0       . 


This  has  the  two  solutions 

\^{z)   =  a(z)  -  i/<(z)\ 
U5)  >><(z)  =/t(z)-a^(z) 

X^Cz)  =  a(z)  +  i/«(z)J 

The  conditions  (a)  and  (h)  of  section  II  are  extended  to  exclude  all  values  of 
z  such  that  "b(z)  =  a  (z).  The  matrix  which  diagonalizes  A(z)  is  found  to  he 


■    -1 


and  time 


-1&- 


1+7) 


yo 


>Vj^(z) 


R(z) 


X^  *  (Aj^(z)-X^(z)) 


rx2(o\(z) 


J-X^UTTI 


Xj^^) 


(^.(zUJO) 

1     2 


In  equation  (Uy)  and  the  following  the  dots  denote  differentiation.   The 
solution  developed  in  section  III  "oecomes 


with  a  sinilar  expression  for  u-C^).  Here  thel)^^(z)  are  expressed  "by- 
equation  (15)  v^here 

m^j^(z)  =  022^^^  ""  ° 


ij' 


U9) 


i2i(r,)  =  X^(.^)ni,2^z) 


™      ^    ^  1  1  r      Miz)  ^   ,    a(z)l 


^<I^^^*Ia%}1^ 


and  iiUp(z)   is  the  complex  conjugate  of  m^Az), 

If   in   equation  (M^  we   taJce  a(z)   =  0   the   resulting   equation, 
7^(2)  +  "bCz)  1/^(2)  =  0,  has  many  important  physical  applications.     In  particular 
it  represents   the  reduction  of  Maxwell's  equations  for  plane  electromagnetic 
vfaves  in  an  isotropic  plane  stratified  inhomogeneous  medium.      The  solution  now 
■becomes 


z 


fM 


Jji^TJ 


[^ 
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n'^V^i2(')-j 


V-^K)) 


50) 


r/«(y)dy 


7/^ 


[n^(z)B,*/:i^^(z)B^]  y^T 


^ 


T-fhere   the   \l.  .(z)   are  defined  as  "before,    with 


51) 


i^^(z)  =  m^gCz)  =  0 


V''>=^[lj^»/<(^>^'j 


i2f  ;t^(y)dy 


■12  f //(y)dy 


mp,U)  =  -//(z   ) 


l.V')^'l 


and/^  (z)  =  J^TzJ     , 

If  the  integrals  in  the  il^,(z)  are  small  then  it  is  clearly  shovm 
that  the  leading  term  of  the  ahove  solution  is  the  ordinary  V/.B.K.  solution. 
Also  the  usual  turning  points  occur  when  U{z)  =   0  which  correspond  to  eoual 
roots  in  the  present  treatment.  Considering  the  electromagnetic  prohlem  men- 
tioned at  the  "beginning  of  this  section,  the  physical  interpretation  hecomes 
quite  clear.   The  positive  and  negative  exponentials  represent  "rising"  and 
"downcoming"  waves  respectively.  A  pc.rticular  Al..(z),  say  -^^it^z),  is 
explicitely 

^   i2jMy)dy.^  ^   -i2?/^(y)ay 

^ii^^)-*£ldi^^^/^K4  ^°  jMj-'-A  °    ' 


2^* 
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This  series  represents  all  contribationB  of  even  numbers  of  continuous  multiple 
reflections  to  the  upgoing  v/ave.  The  term  jfl  npCz)  yields  the  contri'outions  to 
the  u^pgoing  v/ave  of  all  such  odd  numhers  of  continuous  reflections.   Similarly 
for  the  remaining  A  2  ^Az),      Thus  we  can  define  a  continuous  reflection  co- 
efficient p(z)  "by: 

'  1 

52)  /(z)  =  |^/ny4(z)2   , 

The  solution  of  the  present  section  is  very  similar  to  that  obtained  hy  Bremmer 
for  the  electromagnetic  prohlem  [loj. 
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VIII.   Calculation  of  Connection  Matrix. 

The  solution  of  the  previous  section  fail?,  v.-hen  Xt(z)  =  X-Cz).   In 
the  case  of  a(z)  f  0,  this  occurs  for  Td(z)  =  0.   Let  this  turning  point  "be 
located  at  z  =  0  and  in  the  neigh"borhcod  of  this  point  we  take 

53)  ^5(z)  =  k^  z"   ,  «  >  0  . 

Then  the  A  matrix  of  the   equivalent   first   order  system  is 

51+) 


A(z)  =    (  I 

\-kV  0      / 


The  solution  for  z  >  0  and  ^  <  0  is  given  by  equation  (50)  with  the  appropriatt 
B^  in  each  case.   These  solutions  are  connected,  as  explained  in  section  IV, 
"by  the  matrizant 

t(f.z) fUt^u)!   =  s}  {a(x)]  Ji  U(^)1 

55)  ^         f_i 

=  Q     [m.)]    D   {a(x))  .    . 

Here  v;e  have  us^d  the  identities    (g)   and  (9).      Thus  we  need  compiute   only 
z 
Xl        |a(x))   ,   whose   elements   are,    fron  equations   (53)    and    (2) 

z     x^  *     ?U     ^^"5     ^2 


0000 


^3  ^  ^^2^-^^k  -  * 
Q^{zhl-k^}     J    x^dx^dx^+k^J    J      J      Jx^'x^dx^^dxcbcd^i^-* 


0       0  0  0 


56)  -  a     .      ^ 

The  first  two  of  the  above  3eries  yield  UDcn  integration 
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V,2     (<*+2)^ 


11 


(z)  =i:(-ir 


nl{ci*2)      jr[m(oc+2)-l] 


57) 


Applying  the   identities  [ll"\ 


n  n  /  a-1 
a 


-^L        J  p   (a::ij 


n-l  •'^-1P(1* 

nFl  Wl"  ^ 


n) 


to  equations   (j?)  we  o"btain  \llj 

V        -  ~ 

n^j^cz)  =  r(i-v)(vk)  z2j_^  (2vk2^'') 
%(z)  =  r  ( v)(vk)^-^  z^  ^^^  j^-i^'^^^^^ 

From  the  last  tv^o   of  equations   (56  )   and  equations   (58)  there  resiilts 

1(1- ->>    )  1_ 

fl      (2)  =  -kr(l-Y)(Vkr'     z^      ^       J,      ,(2Vkz2^) 


a^''     -^,  x^  .,v .-.,     .  -j^^ 

59) 


ni2^z)  =  ^  P(V)(Vk)^"^     z2     J^     (2>kz2^    ) 


The  elements  of    Jt  L.      \  ■^\  ^^^  easily  computed  "by  finding  the  inverse 
of     LL   f  Aland  replacing  z  ty  f  ,     Then  we  obtain  the  connection  matrix  of 
equation  (55)  as 


1  «.n 

t(r.z)=[ 

m+l  ^  1 

Here  v/e  have  neglected  the  constant  factor  P  (l-V)^  (V)(vk)  and  for  the 
arguments  of  the  Beasel  functions  the  atbreviated  notation 

2^ 
61) 

has  "been  used.   This  matrix  furnishes  the  exact  connection  between  rlgoroTis 
solutions  on  either  side  of  a  turning  point  (of  the  form  given  in  equation  (53)). 
Wiat  are  usually  referred  to  as  connection  formulae  are  concerned  with  connecting 
asyniptotic  solutions  across  a  turning  point.   These  may  he  obtained  from  the 
matrix  (60)  by  using  the  asymptotic  forms  of  the  Bessel  functions  for  large  k 
and  taking  account  of  the  Stokes  phancnenom  which  arises  since  z'  and  r'  have 
different  arguments  |[  I2], 
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Appondix  I.     Matrizant  of  a   "self-Comnuting."  Katrlx. 

To  prove  identity  (6)  wa  consider  a  "self-commating"  matrix  A(x): 
i.e.    a  matrix   such  that 

62)  A(x')  A(x«)  =  A(x")  A(x') 

for  ar'bitrary  x'  and  x"  in  some  interval.   It  is  assumed  that  the  matrizemt 

of  A(x)  exists.  From  the  definition  (2)  this  may  "be  written  as 

Z       Z  X  ^  z. 

X  f  f^  /'tt-l        f  2 


Bow  consider: 

X 


J     A(pd^  =1:  M()d^  A(z).      J       A 

''x  8=1       ^Z.  Z_ 


^H\ 


8-1 


6U) 


I' 


(pd5 


p-i 


A(x) 


A(x) 


h 


-iP-l 


(C)<i^ 


The  first  line  of  (6^4-)  follov;c   from  the  ordinary  fonnala  for  the  deviative  of 
a  pQv/er  of  a  matrix.     The  last   two  results  are  o'btained  "by  applying  the  property 
(6?.),     The  general  term  in  the   sum  (63)  may  now  Ije  evaluated  with  the  aide  of 
the  final  rssult   in  (6U).      Thib  yields 


65) 


4  uw]  =1*  E  ^ll^H^^^l" 

0  n— 1  I      0  * 

X 

=  exp   I       A(pd5 


Appendix  II.  Non-Simple  Elementary  Divisors. 

The  cpce  of  raultiole  latent  roots  may  be  treated  quite  generally  as 

follows.  Let  there  be  r  roots  each  of  multiplicity  s.,  i=l,2, ...,  r;  so  that 
r 

Y~    3.  =  n.  Or  in  other  ^rords,  the  r-elementary  divisors  of  the  matrix^A-XI  ) 
1=1  ^  \     ' 

have  orders  s^^,  S2 ,  s^  respectively.   In  thio  case  it  is  known  that  there 

exists  a  non-singular  matrix  P(z)  such  that 

66)        P"^(0  A(j)  p(0  =.  B(z)  s     ^2.        . 

\  O      3   / 

P 

Here  each  B,(z)  is  of  order  s.  with 


67) 


and  X.(z)  is  a  root  of  the  characteristic  equation  (not  all  X. (z)  need  he 
distinct). 

We  use  the  ahove  matrix  P(z)  in  the  solution  given  hy  equation  (lO) 
and  then  apply  identity  (j)  to  ohtaln 


681) 


U(z)=P(z)  A       {    B+D^  Q        \m(x)]  P'^Cz^)!!^ 


Again  v/e  have  introduced: 

g^j  p-^(z)     ^1^    =     R(z)  S  D(z)  -^NCz)      ,        and 


o  ~ 

M(z)  5       n      \'B+Ii\    H(z)      -Tl       [  »^°^ 


where  D  is  diagonal  and  N  has  zero  diagonal  elements.  From  the  definition  of  the 
ma.trizant  it  is  easily  proven  that: 
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70) 


'4  [Vl\ 


o 


4      ^^2^4 


0 

Here   the  D.  (z)    are   (diagonal)   culs-niatricee   of  D,    each  of   the   same  order  as   the 
associated  B.(z),     The  identity   (7)  applied  to   each  of  the  matrizants  on  the 
right  of  equation  (70)  gives 

0  0  0       ».  o  J 

The  matrizant  of  B.(z),  which  is  given  "by  equation  (67),  can  be  calculated  frcm 
the  definition;  this  yields 


72) 


4  {\^4 


!  \^^)^i 


1 


2(2-2^) 


(z-z^) 


(z-z  ) 

0 


o 


*  1 


Kow  actuations  (70)-(72)  can  he  used  in  the  solution  (68),  "but  the  resulting 
form  of  the  solution  will  contain  essentially  two  matrizants  which  must  he 
evaluated.  It  will  he  recalled  that  in  the  case  of  complete  dlagonalizatioa 
(simple  elementary  divisors)  the  final  form  of  the  solution  involved  the  evalu- 
ation of  only  one  matrizant. 
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Appendljc  III.      Blements  of     P"-^     ^ 
dz 

Let   the   elemente   of  a  non-aingular  matrix  P(z)  "be  p.  ^(z).      Then   the 
elements  P. .(z)   of   the  inverse  matrix  P~   (z)   are 

73^  ^ij^^^  =  iTferr  \i^^^ 


where 


v^'H^aH  •(;;; 


=  cofactor  of  p.j^(z)  in  |p(z)l  , 


The  elements  R,  .(z)  of  the  matrix  (-P~  ^)  are  thus 
ij  dz 


n   cLp„/z) 


«8=1 

Hov/ever,   from  the  definition  (73): 


,      n       dp     (z) 


76) 


rP3j(^^  v=^^=  Ki^'^^^^ 


Uain€  (77 ^  in  (75^  we  oUain 


fe       dz         is  ij     dz  ^^    ^sj        dz 


(78)  V''  =  -  ^ij  ^'»'^<'"  *  ifbr  f,^j  ^  • 

Letting  i  =  j   the   result  of   equations   (22)   and  (23)  follov^ 


^  28  - 
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